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Abstract 

A  simulation  tool  for  modeling  planar  solid  oxide  fuel  cells  is  demonstrated.  The  tool  combines  the  versatility  of  a  commercial 
computational  fluid  dynamics  simulation  code  with  a  validated  electrochemistry  calculation  method.  Its  function  is  to  predict  the  flow  and 
distribution  of  anode  and  cathode  gases,  temperature  and  current  distributions,  and  fuel  utilization.  A  three-dimensional  model  geometry, 
including  internal  manifolds,  was  created  to  simulate  a  generic,  cross-flow  stack  design.  Similar  three-dimensional  geometries  were  created 
for  simulation  of  co-flow,  and  counterflow  stack  designs.  Cyclic  boundary  conditions  were  imposed  at  the  top  and  bottom  of  the  model 
domains,  while  the  lateral  walls  were  assumed  adiabatic.  The  three  cases  show  that,  for  a  given  average  cell  temperature,  similar  fuel 
utilizations  can  result  irrespective  of  the  flow  configuration.  Temperature  distributions  however,  which  largely  determine  thermal  stresses 
during  operation,  are  dependent  on  the  chosen  design  geometry/flow  configuration.  The  co-flow  case  had  the  most  uniform  temperature 
distribution  and  the  smallest  thermal  gradients,  thus  offers  thermo- structural  advantages  over  the  other  flow  cases. 
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1.  Introduction 

The  planar-type  solid  oxide  fuel  cell  (SOFC)  is  desirable 
for  many  applications  due  to  its  compactness.  Cell  active 
areas  can  be  stacked  in  close  proximity  to  yield  large  power 
densities.  Development  of  the  planar  SOFC  has  challenges 
related  to  the  high-temperature  operation  (600-1000  °C). 
Many  designs  use  excess  cathode  airflow  to  remove  heat 
during  steady-state  operation.  The  distribution  and  utilization 
of  fuel  over  the  anode  and  the  distribution  of  relatively  cooler 
cathode  air  interact  to  determine  the  temperature  gradients. 
One  of  the  many  challenges  associated  with  developing  this 
type  of  SOFC  is  to  minimize  the  non-uniform  temperature 
distributions  that  contribute  to  stress  in  the  planar  compo¬ 
nents  (stress  analyses  are  not  reported  in  this  paper). 

Fuel  utilization  and  average  cell  temperature  are  con¬ 
trolled  by  the  delivery  rate  and  the  temperatures  of  fuel  and 
air  to  the  cell.  Fuel  distribution  is  critical  due  to  the 
exothermic  electrochemical  reactions:  cell  areas  covered 
by  concentrated  fuel  become  electrochemically  active,  lead¬ 
ing  to  increased  local  temperatures,  which  yields  even  faster 
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reaction  rates.  Cell  areas  covered  by  depleted  fuel  become 
inactive,  thus  cooler,  which  leads  to  even  slower  localized 
reaction  rates.  Increased  fuel  flow  tends  to  increase  uni¬ 
formity  of  the  reaction  rates  across  the  active  area  and  to 
decrease  fuel  utilization.  Decreased  fuel  flow  tends  to 
increase  fuel  utilization,  but  can  cause  local  fuel  depletion 
and  cold  spots  that  exacerbate  temperature  non-uniformities. 
Air  and  fuel  inlet  temperatures  also  affect  the  reaction  rates, 
cell  temperature,  and  the  fuel  utilization.  Therefore,  man¬ 
agement  of  the  flow  of  air  and  fuel,  and  the  distribution  of 
each,  is  critical  to  stable  operation  of  the  cell.  Other  factors 
that  affect  the  temperature  distribution  and  fuel  utilization 
are  the  relative  orientations  of  the  fuel  and  airflows,  the  mass 
and  dimensions  of  the  cell  components,  and  dimensions  of 
flow  channels.  Thus,  geometry  is  important  in  establishing  a 
well-operating  cell.  In  order  to  efficiently  develop  and 
optimize  planar  SOFC  stacks,  it  is  convenient  to  have  the 
capability  to  experiment  numerically  with  the  effects  of 
geometric  configurations  on  operation  and  performance. 

Modeling  of  the  planar  SOFC  during  steady  operation  to 
calculate  temperatures  and  stresses  is  ongoing  and  becoming 
increasingly  three-dimensional  [1-3].  Investigations  of  pla¬ 
nar  SOFC  operation  [4]  and  planar  SOFC  performance  [5] 
have  predicted  cell  temperature  and  electrical  current  dis¬ 
tributions  for  various  flow  configurations.  The  purpose  of  the 
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Nomenclature 

FhJ  diffusional  energy  flux  in  direction  xj 
yfg  determinant  of  metric  tensor 

h  enthalpy 

/ave  average  electrical  current,  A/cm2 

mm  mass  fraction  of  mixture  component  m 
Mm  molecular  weight  of  mixture  component  m 
p  piezometric  pressure 

R  universal  gas  constant 

sh  enthalpy  source  term 

Si  momentum  source  components 

sm  mass  source  term 

T  absolute  temperature,  Kelvin 

rPEN  average  cell  (or  PEN)  temperature 

A7pEN  difference  between  maximum  and  minimum 
temperature  on  active  area  of  cell 
ut  fluid  velocity  component  in  direction  xt 

Uj  fluid  velocity  relative  to  coordinate  system 

moving  with  velocity  ucj 

V  voltage 

xt  Cartesian  coordinate  (i  =  1,  2,  3) 

Greek  Letters 
p  density,  kg/m3 

t ij  stress  tensor  component 


present  work  is  to  demonstrate  a  modeling  tool  to  predict  the 
fuel  utilization,  electric  current,  and  temperature  distribu¬ 
tions  in  planar  SOFC  stacks  with  arbitrary  three-dimensional 
geometries.  For  this  study,  cross-,  co-,  and  counterflow 
configurations  are  examined. 

The  modeling  tool  couples  a  validated  electrochemistry 
calculation  method  with  a  commercial  computational  fluid 
dynamics  (CFD)  simulation  code.  In  the  calculations,  the 
code  solves  the  finite-volume  Navier-Stokes  and  transport 
equations  to  obtain  the  flow,  species  concentrations,  and 
temperature  at  each  location  in  the  cell/stack.  This  informa¬ 
tion  is  passed  to  the  electrochemistry  module  (subroutine). 
Here  the  local  current  density  is  calculated,  based  on  the 
temperature,  applied  cell  voltage,  and  pressures.  The  current 
is  then  used  to  calculate  the  hydrogen  combustion  rate.  The 
water-gas  shift  reaction  rate  is  assumed  to  maintain  equili¬ 
brium  conditions.  Heat  generation  rates  and  species  source 
(fuel/air)  rates  are  supplied  to  the  code  based  on  the  hydro¬ 
gen  combustion  and  shift  reaction  rates.  Gas  species  con¬ 
centrations  and  temperature  distributions  are  then  calculated 
for  the  next  iteration,  and  so  on,  until  equilibrium  is  achieved 
for  the  shift  reaction. 


2.  Modeling  approach 


Fig.  1.  Schematic  view  of  a  typical  cross-flow  planar  interconnect  plate. 


electrode  (PEN)  and  an  interconnect  plate  “stacked” 
together.  Fig.  1  shows  a  schematic  of  an  interconnect  plate. 
The  schematic  is  displayed  as  being  transparent  to  show  top 
and  bottom  features.  In  this  example,  the  interconnect  is 
formed  to  create  internal  gas  flow  manifolds  and  flow 
passages  above  and  below  the  plate  in  a  cross-flow  config¬ 
uration.  The  active  cell  area  in  Fig.  1  is  that  formed  by  the 
intersection  of  air  and  fuel  flows  in  the  X—Y plane.  A  cell  unit 
(or  single-cell  stack)  is  formed  when  interconnects,  similar 
to  that  in  Fig.  1,  are  placed  above  and  below  a  PEN.  Fig.  2 
shows  a  partial  schematic  view  of  a  vertical  slice  through  a 
Y-Z  plane  of  a  multiple-cell,  cross-flow  stack.  As  in  Fig.  1, 
air  flows  along  the  T-direction.  Fuel  flow  is  along  the 
X-direction.  Fuel  flow  vectors  are  orthogonal  to  the  page 
and  are  represented  in  Fig.  2  by  X  in  the  fuel  region.  The 
location  of  a  cell  unit  along  Z  is  shown  at  left  in  the  figure. 

A  generic  single-cell,  cross-flow,  stack  model,  including 
internal  manifolds,  was  created  for  benchmark  calculations 
of  the  fuel  utilization,  gas,  temperature,  and  electrical 
current  distributions.  The  model  represented  a  repeating 
cell  unit  residing  in  the  center  of  a  large  stack  such  that 
no  end  effects  were  present  at  the  cell  of  interest.  The  mesh 
contained  23  elements  along  the  X-  and  ^directions  (529 
elements)  to  resolve  a  rectangular  1 16.6  cm2  active  cell  area 
on  the  PEN.  The  760  mm  thick  PEN  was  modeled  as  a 
monolithic  structure  containing  one  element  in  the  Z-direc- 
tion.  Fuel  and  air  passages  above  and  below  the  PEN  each 
were  subdivided  into  three  elements.  Individual  flow  chan¬ 
nels  created  by  webs  in  the  interconnect  (as  in  Fig.  1)  were 
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The  repeating  unit  of  a  typical  planar  SOFC  stack  is 
constructed  of  a  positive  electrode-electrolyte-negative 


Fig.  2.  Schematic  section  view  of  a  typical  cross-flow  stack.  X  in  the  fuel 
region  represent  fuel  flow  vectors  perpendicular  to  the  page. 
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not  modeled.  Instead,  we  assumed  that  the  region  was  one 
large  flow  passage  and  a  material  placed  in  the  region  would 
act  as  current  collector  and  make  the  electrical  contact 
between  PEN  and  interconnect.  A  glass  seal  between  the 
PEN  and  interconnect  was  included  in  the  model.  The  three- 
dimensional  model  containing  one  repeating  cell  unit  (with 
PEN,  air  and  fuel  channels,  glass  seal,  and  1/2  interconnect 
thickness  top  and  bottom,  in  the  Z-direction)  including 
internal  gas  flow  manifolds  contained  26,255  computational 
elements.  Mesh  resolution  was  chosen  based  on  simulation 
run-time/accuracy  tests  and  represents  an  optimized  mesh. 

2.7.  Thermo-fluid  model 


The  widely  used  STAR-CD  [6]  code  was  selected  for  the 
thermo-fluids  calculations  and  as  host  for  the  electrochem¬ 
istry  calculations.  The  electrochemistry  calculations  will  be 
discussed  in  a  later  section.  STAR-CD  is  a  CFD  code 
that  solves  the  algebraic  finite- volume  equations.  For  each 
of  the  computational  mesh  elements,  in  the  fluid  regions, 
STAR-CD  solved  the  Navier-Stokes  equations  for  conser¬ 
vation  of  mass 


1  d  9 

and  momentum 
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In  the  solid  regions,  STAR-CD  conserved  specific  energy. 
The  equation  is  similar  to  (3)  without  the  pressure  terms. 
Both  the  air  and  fuel  streams  were  modeled  as  ideal  gas 
mixtures  with  density  given  by: 


P 


P_ 

RT 


(4) 


where  mm  is  the  mass  fraction  of  a  component  with  mole¬ 
cular  weight  Mm.  Specific  heats  of  the  gaseous  mixtures 
were  calculated  as  mass  weighted  averages  of  the  compo¬ 
nent  values.  Gas  viscosities  and  thermal  conductivities  were 
taken  to  be  constant  at  the  average  cell  temperature. 

Heat  transfer  between  the  fluid  and  solid  materials  was 
limited  to  conduction  and  convection.  The  flow  channels  had 
large  aspect  ratios  (length-to-height,  roughly  100:1  on  the 
cathode  side  and  200:1  on  the  anode  side)  such  that  view 
factors  for  radiant  heat  exchange  downstream  were  small. 
Additionally,  the  geometric  form  of  the  material  used  for 
current  collection  would  obstruct  and  further  reduce  the 


Table  1 

Solid  material  properties 


Cell 

Component 

Density 

(kg/m3) 

Thermal  conductivity 
(W/m  K) 

Specific  heat 
(J/kg  K) 

Interconnect 

7700 

13 

0.8 

PEN 

4200 

2 

0.623 

Glass  Seal 

2500 

1.4 

0.6 

contribution  of  thermal  radiation  in  the  downstream  direc¬ 
tion.  The  effect  of  radiation,  therefore,  was  small  relative  to 
the  other  modes  of  heat  transfer  and  was  neglected  in  the 
calculations. 

Solid  material  properties  used  in  the  calculations  are  listed 
in  Table  1.  Specific  heats  were  reduced  by  factor  103  to 
expedite  the  calculations.  This  is  an  acceptable  practice  for 
steady-state  simulations. 

2.2.  Boundary  conditions 

The  single-cell  stack  model  represented  a  repeating  cell 
unit  in  the  center  of  a  large  stack;  cyclic  boundary  conditions 
were  imposed  at  the  top  and  bottom  of  the  model  domain 
(i.e.  Z-direction  in  Figs.  1  and  2).  Walls  at  the  periphery  of 
the  stack  (at  the  X-  and  7-direction  extremes)  were  assumed 
adiabatic.  Constant  temperature,  velocity,  and  gaseous  com¬ 
position  were  imposed  at  the  inflow  boundaries  for  the  fuel 
and  air. 

2.3.  Electrochemical  model 

The  unit  cell  electrochemistry  model  (UCEM)  predicts 
local  electrical  responses  to  changes  in  fuel  flow  rate,  local 
fuel  composition,  and  local  temperature.  The  response 
characteristics  are  adjustable  to  component  materials  and 
dimensions  as  well  as  to  electrode  porosity.  The  model 
accounts  for  overpotentials  due  to  Ohmic  resistance  of 
the  cell  components,  contact  resistance  at  electrode-current 
collector  interfaces,  charge  transfer  at  the  electrodes,  and 
diffusion  of  reactants  into  and  products  out  of  the  porous 
electrodes.  The  UCEM  predicts  the  current-voltage  (7-V) 
response  curve  based  on  the  local  conditions  for  a  specific 
unit  area  within  the  active  area  of  the  cell.  The  model 
contains  adjustable  parameters,  which  are  set  by  calibrating 
the  model  to  I-V  data  obtained  from  operation  of  a  small 
(1  in.  diameter)  single-cell  SOFC  under  a  range  of  fuel 
compositions  and  temperatures.  Construction  and  calibra¬ 
tion  of  the  UCEM  is  described  elsewhere  [7,8]. 

2.4.  Electrochemistry  equations  in  STAR-CD 

The  interface  between  STAR-CD  and  the  electrochemical 
(EC)  model  was  implemented  via  the  POSDAT  routine  in 
STAR-CD.  POSDAT  was  called  after  each  STAR-CD  itera¬ 
tion.  The  EC  model  was  called  from  within  POSDAT  for 
each  PEN  element  in  the  model  domain.  Information  passed 


112 


K.P.  Recknagle  et  al./ Journal  of  Power  Sources  113  (2003)  109-114 


included  temperature,  gas  concentrations  (in  the  adjacent 
fluid  streams),  and  PEN  voltage.  Voltage,  rather  than  current, 
is  held  constant  across  a  PEN  to  facilitate  multiple-cell  stack 
simulations.  A  non-linear  root  solver  was  used  to  find  the 
local  electrical  current  on  the  PEN  corresponding  to  the  PEN 
voltage.  The  local  current  was  found  according  to  the 
UCEM.  The  Nernst  voltage  and  overvoltages  due  to  But¬ 
ler- Volmer,  resistive,  and  concentration  polarizations  are 
also  calculated.  A  thermodynamic  equilibrium  constant  for 
the  water-gas  shift  reaction  is  then  calculated  for  each  PEN 
element  using  Gibb’s  free  energies  and  the  local  tempera¬ 
ture.  This  is  compared  to  the  local  chemical  equilibrium 
constant  based  on  the  instantaneous  concentrations  of  each 
gas  species  in  the  reactions.  Agreement  between  the  two 
equilibrium  constants,  to  within  a  small  error,  denotes 
convergence  of  the  STAR-CD/EC  simulation.  Local  PEN 
heat  generation  is  also  calculated  during  each  iteration  and 
provided  to  STAR-CD  for  the  next  iteration. 

2.5.  Calculations  and  results 

Conditions  for  the  benchmark  cross-flow  case  were  cho¬ 
sen  to  achieve  average  cell  temperature  and  fuel  utilization 
of  roughly  750  °C  and  60-70%,  respectively,  for  a  cell 
voltage  of  0.7  V.  Air  was  delivered  to  the  cathode  at 
0.25  gm/s  and  625  °C.  The  benchmark  fuel  was  a  POx- 
reformed  gasoline  with  45%  anode  recycle,  delivered  to 
the  anode  at  0.018  gm/s  (9.5E  —  4  mol/s)  and  625  °C.  Molar 
composition  of  the  fuel  was  36%  H2,  35%  CO,  5%  H20,  5% 
C02,  and  19%  N2. 

The  resulting  steady- state  cell  (PEN)  average  temperature 
was  769  °C  with  maximum  and  minimum  PEN  temperatures 
of  643  and  912  °C,  respectively.  Fuel  utilization  was  61.7%. 
The  PEN  temperature  distribution  for  this  case  is  presented 
in  Fig.  3.  The  temperature  was  highest  near  the  fuel  inlet  and 
air  outlet.  The  Nernst  potential  and  reaction  rates  on  the  PEN 
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Fig.  3.  Predicted  PEN  temperature  distribution  for  single-cell  stack  in  cross- 
flow  configuration.  Temperature  range:  643-912  °C  (769  °C  average). 


Table  2 


Operating  conditions  for  single-cell  stack  simulations  (cell  voltage  =  0.7  V) 


Flow  configuration 

Air  delivery 

Fuel  delivery 

gm/s 

°C 

gm/s 

°C 

Cross-flow 

0.25 

625 

0.018 

625 

Co-flow 

0.25 

625 

0.018 

625 

Counterflow 

0.25 

595 

0.018 

595 

were  highest  near  the  fuel  inlet  and  decreased  in  the  direc¬ 
tion  of  the  fuel  flow.  Air  flowing  past  the  PEN  was  most 
effective  at  cooling  near  the  air  inlet  and  carried  heat 
generated  in  the  PEN  toward  the  air  outlet,  thus  moving 
the  “hot  spot”  toward  the  air  outlet. 

An  attempt  to  achieve  a  more  uniform  temperature  dis¬ 
tribution  led  to  an  investigation  into  the  effect  of  cell 
geometry  on  the  fuel/air  distributions  and  fuel  utilization. 
In  addition  to  the  cross-flow  case,  cases  with  co-flow  and 
counterflow  configurations  were  generated.  Active  area, 
flow  passage  heights,  and  interconnect  dimensions  from 
the  cross-flow  case  were  all  retained.  Operating  conditions 
for  these  other  cases  were  adjusted  to  roughly  reproduce  the 
fuel  utilization  and  average  cell  temperature  obtained  by  the 
cross-flow  case.  Conditions  (delivery  rates  and  temperatures 
of  the  air  and  fuel  streams)  are  listed  in  Table  2. 

Results  from  the  three  flow  configuration  cases  are  shown 
in  Figs.  4-6,  and  listed  in  Table  3.  Temperatures  are  shown  in 
Fig.  4.  Current  density  distributions  are  shown  in  Fig.  5. 
Mass  concentration  of  the  hydrogen  is  shown  in  Fig.  6.  Cell 
temperature  difference,  average  cell  temperature,  current 
density,  and  fuel  utilization  are  itemized  in  Table  3. 

The  cross-flow  PEN  temperature  distribution  on  the  active 
cell  area  is  shown  in  plan  view  in  Fig.  4a.  Also  in  Fig.  4  are 
the  cell  area  temperatures  for  the  co-flow  (Fig.  4b)  and 
counterflow  (Fig.  4c)  cases.  In  each  case  the  temperature 
increased  along  the  airflow  direction  (top  to  bottom  in 
Figs.  4a-c),  reaching  a  maximum  near  the  air  exit.  Of  the 
three  flow  configurations  tested,  the  co-flow  case  had  the 
most  uniform  temperature  distribution  and  smallest  tem¬ 
perature  difference  (AT)  from  air  inlet  to  outlet.  This  is  due 
to  offsetting  effects  of  air  near  the  inlet,  at  its  coolest,  being 
aligned  with  the  fuel  inlet  where  the  Nernst  potential  is 
highest.  Cross-  and  counterflow  cases  had  ATs  of  similar 
magnitude  across  the  PEN  (see  Table  3),  both  larger  than  that 
of  the  co-flow  case  by  more  than  80  °C. 


Table  3 

Single-cell  stack  simulation  results  (cell  voltage  =  0.7  V) 


Flow 

configuration 

ATpen 

(°C) 

TpEN 

(°C) 

Aive 

(A/cm2) 

Fuel 

utilization  (%) 

Cross-flow 

269 

769 

0.69 

62 

Co-flow 

184 

763 

0.71 

64 

Counterflow 

267 

758 

0.73 

63 
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(a)  Cross-Flow  (b)  Co-Flow 


Fig.  4.  Temperature  distribution  on  the  active  cell  area,  °C. 


Consistent  with  the  temperature  distribution,  the  current 
density  (see  Fig.  5)  was  most  uniform  in  the  co-flow  case.  In 
the  cross-flow  and  counterflow  cases,  current  density  was 
high  near  the  fuel  inlet  and  air  outlet  corresponding  to  high 
temperature  and  Nernst  potential  in  that  region.  This  result 
was  not  present  in  the  co-flow  case  where  the  maximum 
current  density  was  roughly  midstream  from  inlet  to  exit — 
again  due  to  the  coincident  location  of  air  and  fuel  inlets  that 
is  a  characteristic  of  the  co-flow  configuration. 

Mass  concentration  of  H2  is  shown  for  the  three  cases  in 
Fig.  6.  The  primary  difference  was  the  two-dimensional 
nature  of  the  cross-flow  case.  The  co-flow  and  counterflow 
cases  had  low  fuel  concentration  across  the  full  width  of  the 


fuel  outlet  creating  a  “line”  minimum.  The  cross-flow  case 
had  a  “point”  minimum  that  was  lower  in  concentration 
than  either  of  the  other  two  cases.  When  trying  to  maximize 
fuel  utilization,  a  “line”  minimum  (rather  than  a  “point” 
minimum)  would  be  more  suitable  to  decrease  the  risk  of 
localized  fuel  depletion.  Preventing  fuel  depletion,  and 
subsequent  cold  spots,  contributes  to  smaller  thermal  gra¬ 
dients  and  less  severe  thermal  stresses  thus  helping  to  ensure 
the  structural  integrity  of  the  cell  components. 

Table  3  lists  the  cell  (or  PEN)  AT,  the  average  PEN 
temperature,  current  density,  and  fuel  utilization  for  the 
three  flow  configuration  cases.  Of  interest  is  the  roughly 
30%  reduction  of  the  AT  in  the  co-flow  case  compared  to  the 


Fig.  6.  Hydrogen  distribution  in  the  fuel  stream  adjacent  to  the  active  cell  area.  Mass  percent  H2,  kg/kg  x  100%. 
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other  flow  cases.  This  occurs  while  average  PEN  tempera¬ 
ture,  current  density,  and  fuel  utilization  differ  by  no  more 
than  3%  between  the  co-flow  and  either  of  the  other  two 
cases.  The  co-flow  case,  due  to  its  relatively  uniform  tem¬ 
perature  distribution  and  smaller  thermal  gradients,  offers 
thermo-structural  advantages  over  the  other  flow  cases. 

3.  Summary  and  conclusions 

A  simulation  tool  for  three-dimensional  modeling  of  planar 
solid  oxide  fuel  cells  has  been  demonstrated.  The  tool  com¬ 
bines  a  commercial  CFD  simulation  code  and  a  validated 
electrochemistry  calculation  method,  enabling  the  calcula¬ 
tions  to  be  performed  on  complex,  three-dimensional  meshes. 
The  effects  of  cell  flow  configuration  on  the  distribution  of 
temperature,  current  density,  and  fuel  species  was  investi¬ 
gated.  Three  flow  configurations,  cross-flow,  co-flow,  and 
counterflow  were  examined  using  identical  inflow  rates  and 
variable  inflow  temperatures.  It  was  found  that,  for  similar 
fuel  utilization  and  average  cell  temperature,  the  co-flow  case 
had  the  most  uniform  temperature  distribution  and  the  smal¬ 
lest  thermal  gradients,  which  is  advantageous  with  respect  to 
ensuring  structural  integrity  of  the  stack  components. 
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